Software Vault: The Diamond Collection

home *** CD-ROM | disk | FTP | other *** search

/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / NEWMATA.TXT < prev next >

Wrap

Text File | 1995-01-20 | 71KB | 1,881 lines

DOCUMENTATION FOR NEWMAT08, A MATRIX LIBRARY IN C++ Copyright (C) 1991,2,3,4,5: R B Davies This is the "how to use" documentation for newmat. Background information on the structure of the library is given in newmatb.txt. This document is available as an ascii file, newmata.txt, and in hypertext format for reading with "Mosaic". Cross-references in the ascii version are given as section numbers. Introduction 1 Getting started 2 Reference manual 3 Error messages 4 Bugs 5 Files in newmat08 6 Problem report form 7 1 Introduction = ============ Conditions of use 1.1 Description 1.2 Is this the library for you? 1.3 Other matrix libraries 1.4 Where to find this library 1.5 How to contact the author 1.6 Change history 1.7 References 1.8 1.1 Conditions of use === ========== == === Copyright (C) 1991,2,3,4,5: R B Davies. Permission is granted to use and distribute but not to sell except for costs of distribution. Distribution as part of low cost CD-ROM collections is welcomed. ------------------------------------------------------------------------------ Please understand that there may still be bugs and errors. Use at your own risk. I take no responsibility for any errors or omissions in this package or for any misfortune that may befall you or others as a result of its use. ------------------------------------------------------------------------------ Please report bugs to me at robert.davies@vuw.ac.nz or Compuserve 72777,656 When reporting a bug please tell me which C++ compiler you are using (if known), and what version. Also give me details of your computer (if known). Tell me where you downloaded your version of my package from and its version number (eg newmat03 or newmat04). (There may be very minor differences between versions at different sites). Note any changes you have made to my code. If at all possible give me a piece of code illustrating the bug. See the problem report form [7]. "Please do report bugs to me." 1.2 General description === ======= =========== The package is intended for scientists and engineers who need to manipulate a variety of types of matrices using standard matrix operations. Emphasis is on the kind of operations needed in statistical calculations such as least squares, linear equation solve and eigenvalues. It supports matrix types Matrix (rectangular matrix) nricMatrix (variant of rectangular matrix) UpperTriangularMatrix LowerTriangularMatrix DiagonalMatrix SymmetricMatrix BandMatrix UpperBandMatrix (upper triangular band matrix) LowerBandMatrix (lower triangular band matrix) SymmetricBandMatrix RowVector (derived from Matrix) ColumnVector (derived from Matrix). Only one element type (float or double) is supported. The package includes the operations *, +, -, concatenation, inverse, transpose, conversion between types, submatrix, determinant, Cholesky decomposition, QR triangularisation, singular value decomposition, eigenvalues of a symmetric matrix, sorting, fast Fourier transform, printing and an interface with "Numerical Recipes in C". It is intended for matrices in the range 15 x 15 to the maximum size your machine will accomodate in a single array. For example 90 x 90 (125 x 125 for triangular matrices) in machines that have 8192 doubles as the maximum size of an array. The number of elements in an array cannot exceed the maximum size of an "int". The package will work for very small matrices but becomes rather inefficient. A two-stage approach to evaluating matrix expressions is used to improve efficiency and reduce use of temporary storage. The package is designed for version 2 or 3 of C++. It works with Borland (3.1 & 4.0), Microsoft (7 & 8) and Watcom (10) C++ on a PC and AT&T C++ (2.1 & 3) and Gnu G++ (2.3.3, 2.5.8 & 2.6.0). It works with some problems with Zortech C++ (version 3.04). 1.3 Is this the library for you? === == ==== === ======= === ==== Do you 1: understand * to mean matrix multiply and not element by element multiply 2: need matrix operators such as * and + defined as operators so you can write things like X = A * (B + C); 3: need a variety of types of matrices (but not sparse) 4: need only one element type (float or double) 5: work with matrices in the range 10 x 10 up to what can be stored in one block in memory 6: tolerate a large package Then maybe this is the right package for you. 1.4 Other matrix libraries === ===== ====== ========= For commercial packages that support multiple element types, look at M++ from Dyad [phone in the USA (800)366-1573, (206)637-9427, fax (206)637-9428], or the Rogue Wave matrix package [(800)487-3217, (503)754-3010, fax (503)757-6650]. For details of other matrix packages look at the listing from Keith Briggs. A recent version of this is in file kbriggs.txt included with this package. Also look at Ajay Shah's list of free numerical software in C and C++. The file is numcomp-free-c.gz in pub/C-numanal on usc.edu. Nikki Locke produces a list of C++ libraries in general - see pub/usenet-by-group/comp.lang.c++/c++_FAQ/libraries in rtfm.mit.edu. Also look on Compuserve in the forums of the various C++ compilers. 1.5 Where to find this library === ===== == ==== ==== ======= * Compuserve (Borland library, zip format) * SimTel (msdos.cpluspls directory, zip format) - see oak.oakland.edu * comp.sources.misc on Internet (shar format) - see wuarchive.wustl.edu 1.6 How to contact the author === === == ======= === ====== Robert Davies 16 Gloucester Street Wilton Wellington New Zealand "email:" robert.davies@vuw.ac.nz 1.7 Change history === ====== ======= Newmat08 - January, 1995: Corrections to improve compatibility with Gnu, Watcom. Concatenation of matrices. Elementwise products. Option to use compilers supporting exceptions. Correction to exception module to allow global declarations of matrices. Fix problem with inverse of symmetric matrices. Fix divide-by-zero problem in SVD. Include new QR routines. Sum function. Non-linear optimisation. GenericMatrices. Newmat07 - January, 1993 Minor corrections to improve compatibility with Zortech, Microsoft and Gnu. Correction to exception module. Additional FFT functions. Some minor increases in efficiency. Submatrices can now be used on RHS of =. Option for allowing C type subscripts. Method for loading short lists of numbers. Newmat06 - December 1992: Added band matrices; 'real' changed to 'Real' (to avoid potential conflict in complex class); Inject doesn't check for no loss of information; fixes for AT&T C++ version 3.0; real(A) becomes A.AsScalar(); CopyToMatrix becomes AsMatrix, etc; .c() is no longer required (to be deleted in next version); option for version 2.1 or later. Suffix for include files changed to .h; BOOL changed to Boolean (BOOL doesn't work in g++ v 2.0); modifications to allow for compilers that destroy temporaries very quickly; (Gnu users - see the section of compiler performance). Added CleanUp, LinearEquationSolver, primitive version of exceptions. Newmat05 - June 1992: For private release only Newmat04 - December 1991: Fix problem with G++1.40, some extra documentation Newmat03 - November 1991: Col and Cols become Column and Columns. Added Sort, SVD, Jacobi, Eigenvalues, FFT, real conversion of 1x1 matrix, "Numerical Recipes in C" interface, output operations, various scalar functions. Improved return from functions. Reorganised setting options in "include.hxx". Newmat02 - July 1991: Version with matrix row/column operations and numerous additional functions. Matrix - October 1990: Early version of package. 1.8 References === ========== The matrix inverse routine and the sort routines are adapted from "Numerical Recipes in C" by Press, Flannery, Teukolsky, Vetterling, published by the Cambridge University Press. Many of the advanced matrix routines are adapted from routines in "Handbook for Automatic Computation, Vol II, Linear Algebra" by Wilkinson and Reinsch, published by Springer Verlag. 2 Getting started = ======= ======= Overview 2.1 Customising 2.2 Compiler performance 2.3 Updating from previous versions 2.4 Example 2.5 Testing 2.6 2.1 Overview === ======== I use .h as the suffix of definition files and .cpp as the suffix of C++ source files. You will need to compile all the *.cpp files except example.cpp, the tmt*.cpp files, sl_ex.cpp and ex_nl.cpp to get the complete package. Ideally you should store the resulting object files as a library. The tmt*.cpp files are used for testing [2.6], example.cpp is an example [2.5] and sl_ex.cpp and ex_nl.cpp are examples of the non-linear [3.26] solve and optimisation routines. I include a number of "make" files for compiling the example and the test package. See the files section [6] for details. But with Borland and Watcom, its pretty quick just to load all the files in the interactive environment by pointing and clicking. My Borland "make" files are generated from the interactive environment. Use the large or flat model when you are using a PC. Do not "outline" inline functions. You may need to increase the stack size. The "make" files for the Unix compilers link a .cxx file to each .cpp file since some of these compilers do not recognise .cpp as a legitimate extension for a C++ file. I suggest you delete this part of the "make" file and rename the .cpp files to something your compiler recognises. My "make" files for Unix systems are for use with 'gmake' rather than 'make'. Ordinary 'make' works with them on the Sun but not the Silicon Graphics or HP machines. Your source files that access the newmat you will need to #include one or more of the following files. include.h: if you want to access just the compiler options newmat.h: to access just the main matrix library (includes include.h) newmatap.h: to access the advanced matrix routines such as Cholesky decomposition, QR triangularisation etc (includes newmat.h) newmatio.h: to access the output routines (includes newmat.h) You can use this only with compilers that support the standard input/output routines including manipulators. It cannot be used with Zortech or early versions of Gnu. newmatnl.h: to access the non-linear optimisation routines (includes newmat.h) See the section on customising [2.2] to see how to edit include.h for your environment and the section on compilers [2.3] for any special problems with the compiler you are using. 2.2 Customising === =========== The file include.h sets a variety of options including several compiler dependent options. You may need to edit include.h to get the options you require. If you are using a compiler different from one I have worked with you may have to set up a new section in include.h appropriate for your compiler. Borland, Turbo, Gnu, Microsoft, Watcom and Zortech are recognised automatically. If none of these are recognised a default set of options is used. These are fine for AT&T, HPUX and Sun C++. If you using a compiler I don't know about you may have to write a new set of options. Activate the appropriate statement to make the element type float or double. If you are using the standard style for deleting arrays (AT&T version 2.1 or later) make sure Version21 is #defined. If you are using the old style, make sure it is not #defined. There is an option in include.h for selecting whether you use compiler supported exceptions, simulated exceptions, or disable exceptions. Use the option for compiler supported exceptions if and only if you have set the option on your compiler to recognise exceptions. Disabling exceptions sometimes helps with compilers that are incompatible with my exception simulation scheme. I suggest you leave the option TEMPS_DESTROYED_QUICKLY activated, even though the Gnu compiler is the only one I know about that requires it. This stores the "trees" dsecribing matrix expressions on the heap rather than the stack and, surprisingly, seems to give better performance. See the discussion in newmatb.txt for more explanation. Leave the option TEMPS_DESTROYED_QUICKLY_R not activated unless you are using the Gnu G++ [2.3.3] compiler. This option controls whether the ReturnMatrix [3.13] construct uses the stack or the heap. The heap version is rather kludgy and probably should be avoided where possible. The option DO_FREE_CHECK is used for tracking memory leaks and normally should not be activated. Activate SETUP_C_SUBSCRIPTS if you want to use traditional C style element access [3.2]. 2.3 Compiler performance === ======== =========== I have tested this library on a number of compilers. Here are the levels of success and any special considerations. In most cases I have chosen code that works under all the compilers I have access to, but I have had to include some specific work-arounds for some compilers. For the MsDos versions, I use a 486dx computer running MsDos 6 or windows NT. The unix versions are on a Sun Sparc station or a Silicon Graphics or a HP unix workstation. Thanks to Victoria University and Industrial Research Ltd for access to the Unix machines. I have set up a block of code for each of the compilers in include.h. Turbo, Borland, Gnu, Zortech, Microsoft and Watcom are recognised automatically. There is a default option that works for AT&T, Sun C++ 4.0.1 and HPUX. So you don't have to make any changes for these compilers. Otherwise you may have to build your own set of options in include.h. AT&T 2.3.1 Borland 2.3.2 Gnu G++ 2.3.3 HPUX 2.3.4 Microsoft 2.3.5 Sun 2.3.6 Watcom 2.3.7 Zortech 2.3.8 2.3.1 AT&T ===== ==== AT&T C++ 2.1;3.0.1 on a Sun: It works fine with 3.0.1. I haven't been able to test the latest version of the library on 2.1. The previous version worked fine. Except aggregates are not supported in 2.1 and setjmp.h generated a warning message. If you are using "interviews" you may get a conflict with Catch. Either #undefine Catch or replace Catch with CATCH throughout my package. In AT&T 2.1 you may get an error when you use an expression for the single argument when constructing a Vector or DiagonalMatrix or one of the Triangular Matrices. You need to evaluate the expression separately. You cannot run my non-linear [3.26] package with these compilers, because of my use of labels. 2.3.2 Borland ===== ======= Borland C++ 3.1, 4.5: Recently this has been my main development platform, so naturally everything works with this compiler. There was a problem with the library utility in version 2.0 which is now fixed. You will need to use the large model. If you are not debugging, turn off the options that collect debugging information. Make sure you don't run Borland's exceptions and my simulated exceptions at the same time. When running my test program with Borland 4.5 under ms-dos you may run out of memory. Either compile the test routine to run under "easywin" or use simulated exceptions rather than the built in exceptions. Under "easywin" the test program indicates a memory leak. I presume this is partly because of the way windows organises its heap rather than there being a real problem. However there does seem to be genuine memory problem when you use the built-in exceptions. Under "easywin" the automatic clean-up of objects by the exception mechanism does not seem to work. Use my simulated exceptions if this is a problem. 2.3.3 Gnu G++ ===== === === Gnu G++ 2.3.3, 2.5.8: These mostly work. You must enable the options TEMPS_DESTROYED_QUICKLY and TEMPS_DESTROYED_QUICKLY_R. You can't use expressions like Matrix(X*Y) in the middle of an expression and (Matrix)(X*Y) is unreliable. If you write a function returning a matrix, you MUST use the ReturnMatrix [3.13] method described in this documentation. This is because g++ destroys temporaries occuring in an expression too soon for the two stage way of evaluating expressions that newmat uses. Gnu seems to leave some rubbish on the stack. Possibly this is a printer buffer so it may not be a bug. You will have problems with versions of Gnu earlier than 2.3.1. Linux: Gnu G++ 2.5.8 seems a little more touchy than regular G++. But basically the package works. Gnu G++ 2.6.0: Seems OK. There should be better compatibility because version 2.6 retains temporaries until the end of the statement instead of destroying them immediately following the first access. I suggest you enable option TEMPS_DESTROYED_QUICKLY but not TEMPS_DESTROYED_QUICKLY_R. 2.3.4 HP-UX ===== ===== HP 9000 series HP-UX. I have tried the library on two versions of HP-UX. (I don't know the version numbers, the older is a clone of AT&T 3, the newer is HP's version with exceptions). Both worked after the modifications described in this section except that I couldn't compile the non-linear [3.26] package due to my use of labels. With the older version of the compiler I needed to edit the math.h library file to remove a duplicate definition of abs. With the newer version you can set the +eh option to enable exceptions and activate the UseExceptions option in include.h. If you are using my make file, you will need to replace CC with CC +eh where ever CC occurs. I recommend that you do not do this and either disable exceptions or use my simulated exceptions. I get core dumps when I use the built-in exceptions and suspect they are not sufficiently debugged as yet. If you are using my simulated exceptions you may get a mass of error messages from the linker about __EH_JMPBUF_TEMP. In this case get file setjmp.h (in directory /usr/include/CC ?) and put extern in front of the line jmp_buf * __EH_JMPBUF_TEMP; The file setjmp.h is accessed in my file myexcept.h. You may want to change the #include statement to access your edited copy of setjmp.h. 2.3.5 Microsoft ===== ========= Microsoft C++ (7.0, 8.0): Seems to work OK. You must #define TEMPS_DESTROYED_QUICKLY owing to a bug in version 7 (at least) of MSC. There are some notes in the file include.h on changes to run under version 7. I haven't tried the latest version of newmat08 on version 7. Haven't tried visual C++ yet. 2.3.6 Sun ===== === Sun C++ (version 4.0.1): This generally works fine. However, I suspect that there is a problem with my non-linear [3.26] package, even though the program appears to run correctly, probably because of my use of labels. When I set DO_FREE_CHECK, I detect non-existent objects being deleted and the program fails if I use simulated exceptions. 2.3.7 Watcom ===== ====== Watcom C++ (version 10): basically this works fine. Don't try to run Watcom's exceptions and my simulated exceptions at the same time. There does seem to be a problem with expressions such as X = Matrix(A+B)+C; Occasionally the temporary object created by Matrix(A+B) is destroyed twice. In most cases this does not seem to cause a problem. However, I think one should avoid code such as this - in fact, there is generally not much point to such code. Alternatively use my simulated exceptions as the problem seems to occur only with the built in exceptions enabled. 2.3.8 Zortech ===== ======= Zortech C++ 3.04: "const" doesn't work correctly with this compiler, so the package skips all of the statements Zortech can't handle. Zortech leaves rubbish on the heap. I don't know whether this is my programming error or a Zortech error or additional printer buffers. Deactivate the option for version 2.1 in include.h. Does not support IO manipulators. Otherwise the package mostly works, but not completely. Best don't #define TEMPS_DESTROYED_QUICKLY. Exceptions and the nric interface don't work. I think some of the problems are because Zortech doesn't handle conversions correctly, particularly automatic conversions. But, also newmat is just too big for my version of Zortech. Zortech runs much more slowly than Borland and Microsoft. Use the large model and compile as much as possible with optimisation on to save space. You won't be able to get the whole test program to work. Zortech doesn't have io-manipulators so you can't compile the example. I haven't tried the Symantec successors to Zortech. 2.4 Updating newmat08 === ======== ======== Updating from previous 08 betas 2.4.1 Updating from previous versions 2.4.2 2.4.1 Updating from previous 08 betas ===== ======== ==== ======== == ===== If you were using the non-linear optimisation routine in an previous beta, note that rowvectors and columnvectors have been swapped in this version of newmatnl. You now derive from prototype function classes rather than the solution and least squares classes. Simple examples are included with the library. This version includes concatenation operators, elementwise products for matrices and GenericMatrices. See the sections on binary operators [3.6] and unspecified types [3.16]. There is now no need to explicitly set the AT&T option in include.h. I have reorganised the make files for AT&T, Gnu, Microsoft and Watcom. See the files [6] section. 2.4.2 Updating form previous versions ===== ======== ==== ======== ======== This is a minor upgrade on newmat07 to correct errors (one serious) and improve compatibility with various compilers. You should upgrade. * .cxx files are now .cpp files. Some versions of won't accept .cpp. The "make" files for Gnu and AT&T link the .cpp files to .cxx files before compilation and delete the links after compilation. * An option [2.2] in include.h allows you to use compiler supported exceptions, simulated exceptions or disable exceptions. Edit the file include.h to select one of these three options. Don't simulate exceptions if you have set your compiler's option to implement exceptions. * New QR decomposition [3.18] functions. * A non-linear least squares [3.26] class. * No need to explicitly set the AT&T option in include.h. * Concatenation and elementwise multiplication [3.6]. * A new GenericMatrix [3.16] class. * Sum [3.8] function. * Some of the make [6] files reorganised. If you are upgrading from newmat06 note the following: * If you are using << to load a Real into a submatrix change this to =. If you are upgrading from newmat03 or newmat04 note the following * .hxx files are now .h files * real changed to Real * BOOL changed to Boolean * CopyToMatrix changed to AsMatrix, etc * real(A) changed to A.AsScalar() The current version is quite a bit longer that newmat04, so if you are almost out of space with newmat04, don't throw newmat04 away until you have checked your program will work under this version. See the change history [1.7] for other changes. 2.5 Example === ======= An example is given in example.cpp. This gives a simple linear regression example using four different algorithms. The correct output is given in example.txt. The program carries out a rough check that no memory is left allocated on the heap when it terminates. See the section on testing [2.6] for a comment on the reliability of this check and the use of the DO_FREE_CHECK option. I include a variety of make files. Use a command like gmake example -f gnu.mak (Gnu G++) gmake example -f cc.mak (AT&T, HPUX, Sun) nmake example -f ms_nt.mak (Microsoft C++ 8.0) make -f ex_b.mak (Borland C++ 3.1) make -f ex_bc45.mak (Borland C++ 4.5) wmake example.exe -f watcom.mak (Watcom C++ 10A) wmake example.exe -f watco_nt.mak (Watcom C++ 10A) The Borland files were derived from the project file and you will have to edit these to suit your environment. The file ex_bc45.mak is for making a windows NT executable program. The Microsoft file is for the version of C++ that came with the Windows NT 3.1 development kit. The second Watcom file is for making a Windows NT executable. ------------------------------------------------------------------------------ This example uses io manipulators. It will not work with a compiler that does not support the standard io manipulators. ------------------------------------------------------------------------------ 2.6 Testing === ======= The library package contains a comprehensive test program in the form of a series of files with names of the form tmt?.cxx. The files consist of a large number of matrix formulae all of which evaluate to zero (except the first one which is used to check that we are detecting non-zero matrices). The printout should state that it has found just one non-zero matrix. Various versions of the make file (extension .mak) are included with the package. See the section on files [6]. The program also allocates and deletes a large block and small block of memory before it starts the main testing and then at the end of the test. It then checks that the blocks of memory were allocated in the same place. If not then one suspects that there has been a memory leak. i.e. a piece of memory has been allocated and not deleted. This is not completely foolproof. Programs may allocate extra print buffers while the program is running. I have tried to overcome this by doing a print before I allocate the first memory block. Programs may allocate memory for different sized items in different places, or might not allocate items consecutively. Or they might mix the items with memory blocks from other programs. Nevertheless, I seem to get consistent answers from many of the compilers I am working with, so I think this is a worthwhile test. If the DO_FREE_CHECK [2.2] option in include.h is activated, the program checks that each 'new' is balanced with exactly one 'delete'. This provides a more definitive test of no memory leaks. There are additional statements in myexcept.cpp which can be activated to print out details of the memory being allocated and released. Several of the files have a line defining 'REPORT' that can be activated (deactivate the dummy version). This gives a printout of the number of times each of the 'REPORT' statements in the file is accessed. Use a grep with line numbers to locate the lines on which 'REPORT' occurs and compare these with the lines that the printout shows were actually accessed. 3 Reference manual = ========= ====== Constructors 3.1 Accessing elements 3.2 Assignment and copying 3.3 Entering values 3.4 Unary operations 3.5 Binary operations 3.6 Matrix and scalar ops 3.7 Scalar functions 3.8 Submatrices 3.9 Change dimension 3.10 Change type 3.11 Multiple matrix solve 3.12 Memory management 3.13 Efficiency 3.14 Output 3.15 Accessing unspecified type 3.16 Cholesky decomposition 3.17 QR decomposition 3.18 Singular value decomposition 3.19 Eigenvalue decomposition 3.20 Sorting 3.21 Fast Fourier transform 3.22 Numerical recipes in C 3.23 Exceptions 3.24 Cleanup following exception 3.25 Non-linear applications 3.26 3.1 Constructors === ============ To construct an m x n matrix, 'A', (m and n are integers) use Matrix A(m,n); The UpperTriangularMatrix, LowerTriangularMatrix, SymmetricMatrix and DiagonalMatrix types are square. To construct an n x n matrix use, for example UpperTriangularMatrix UT(n); LowerTriangularMatrix LT(n); SymmetricMatrix S(n); DiagonalMatrix D(n); Band matrices need to include bandwidth information in their constructors. BandMatrix BM(n, lower, upper); UpperBandMatrix UB(n, upper); LowerBandMatrix LB(n, lower); SymmetrixBandMatrix SB(n, lower); The integers upper and lower are the number of non-zero diagonals above and below the diagonal (excluding the diagonal) respectively. The RowVector and ColumnVector types take just one argument in their constructors: RowVector RV(n); ColumnVector CV(n); You can also construct vectors and matrices without specifying the dimension. For example Matrix A; In this case the dimension must be set by an assignment statement [3.3] or a re-dimension statement [3.10]. You can also use a constructor to set a matrix equal to another matrix or matrix expression. Matrix A = UT; Matrix A = UT * LT; Only conversions that don't lose information are supported - eg you cannot convert an upper triangular matrix into a diagonal matrix using =. 3.2 Accessing elements === ========= ======== Elements are accessed by expressions of the form 'A(i,j)' where i and j run from 1 to the appropriate dimension. Access elements of vectors with just one argument. Diagonal matrices can accept one or two subscripts. This is different from the earliest version of the package in which the subscripts ran from 0 to one less than the appropriate dimension. Use 'A.element(i,j)' if you want this earlier convention. 'A(i,j)' and 'A.element(i,j)' can appear on either side of an = sign. If you activate the '#define SETUP_C_SUBSCRIPTS' in 'include.h' you can also access elements using the tradition C style notation. That is 'A[i][j]' for matrices (except diagonal) and 'V[i]' for vectors and diagonal matrices. The subscripts start at zero (ie like element) and there is no range checking. Because of the possibility of confusing 'V(i)' and 'V[i]', I suggest you do not activate this option unless you really want to use it. 3.3 Assignment and copying === ========== === ======= The operator = is used for copying matrices, converting matrices, or evaluating expressions. For example A = B; A = L; A = L * U; Only conversions that don't lose information are supported. The dimensions of the matrix on the left hand side are adjusted to those of the matrix or expression on the right hand side. Elements on the right hand side which are not present on the left hand side are set to zero. The operator << can be used in place of = where it is permissible for information to be lost. For example SymmetricMatrix S; Matrix A; ...... S << A.t() * A; is acceptable whereas S = A.t() * A; // error will cause a runtime error since the package does not (yet?) recognise 'A.t()*A' as symmetric. Note that you can not use << with constructors. For example SymmetricMatrix S << A.t() * A; // error does not work. Also note that << cannot be used to load values from a full matrix into a band matrix, since it will be unable to determine the bandwidth of the band matrix. A third copy routine is used in a similar role to =. Use A.Inject(D); to copy the elements of 'D' to the corresponding elements of 'A' but leave the elements of 'A' unchanged if there is no corresponding element of 'D' (the = operator would set them to 0). This is useful, for example, for setting the diagonal elements of a matrix without disturbing the rest of the matrix. Unlike = and <<, Inject does not reset the dimensions of 'A', which must match those of 'D'. Inject does not test for no loss of information. You cannot replace 'D' by a matrix expression. The effect of 'Inject(D)' depends on the type of 'D'. If 'D' is an expression it might not be obvious to the user what type it would have. So I thought it best to disallow expressions. Inject can be used for loading values from a regular matrix into a band matrix. (Don't forget to zero any elements of the left hand side that will not be set by the loading operation). Both << and Inject can be used with submatrix expressions on the left hand side. See the section on submatrices [3.9]. To set the elements of a matrix to a scalar use operator = Real r; Matrix A(m,n); ...... Matrix A(m,n); A = r; 3.4 Entering values === ======== ====== You can load the elements of a matrix from an array: Matrix A(3,2); Real a[] = { 11,12,21,22,31,33 }; A << a; This construction does not check that the numbers of elements match correctly. This version of << can be used with submatrices on the left hand side. It is not defined for band matrices. Alternatively you can enter short lists using a sequence of numbers separated by << . Matrix A(3,2); A << 11 << 12 << 21 << 22 << 31 << 32; This does check for the correct total number of entries, although the message for there being insufficient numbers in the list may be delayed until the end of the block or the next use of this construction. This does not work for band matrices or submatrices, or for long lists. Also try to restrict its use to numbers. You can include expressions, but these must not call a function which includes the same construction. 3.5 Unary operators === ===== ========= The package supports unary operations X = -A // change sign of elements X = A.t() // transpose X = A.i() // inverse (of square matrix A) 3.6 Binary operators === ====== ========= The package supports binary operations X = A + B // matrix addition X = A - B // matrix subtraction X = A * B // matrix multiplication X = A.i() * B // equation solve (square matrix A) X = A | B // concatenate horizontally (concatenate the rows) X = A & B // concatenate vertically (concatenate the columns) X = SP(A, B) // elementwise product of A and B (Schur product) Notes: * If you are doing repeated multiplication. For example 'A*B*C', use brackets to force the order of evaluation to minimize the number of operations. If 'C' is a column vector and 'A' is not a vector, then it will usually reduce the number of operations to use 'A*(B*C)'. * In the equation solve example case the inverse is not explicitly calculated. An LU decomposition of 'A' is performed and this is applied to 'B'. This is more efficient than calculating the inverse and then multiplying. See also multiple matrix solving [3.12]. * The package does not (yet?) recognise 'B*A.i()' as an equation solve and the inverse of 'A' would be calculated. It is probably better to use '(A.t().i()*B.t()).t()'. * Horizontal or vertical concatenation returns a result of type Matrix, RowVector or ColumnVector. * If 'A' is m x p, 'B' is m x q, then 'A | B' is m x (p+q) with the k-th row being the elements of the k-th row of 'A' followed by the elements of the k-th row of 'B'. * If 'A' is p x n, 'B' is q x n, then 'A & B' is (p+q) x n with the k-th column being the elements of the k-th column of 'A' followed by the elements of the k-th column of 'B'. * For complicated concatenations of matrices, consider instead using submatrices [3.9]. 3.7 Matrix and scalar === ====== === ====== The following expression multiplies the elements of a matrix 'A' by a scalar f: 'A * f;' Likewise one can divide the elements of a matrix 'A' by a scalar f:' A / f;' The expressions 'A + f' and 'A - f' add or subtract a rectangular matrix of the same dimension as 'A' with elements equal to 'f' to or from the matrix 'A'. In each case the matrix must be the first term in the expression. Expressions such 'f + A' or 'f * A' are not recognised (yet?). 3.8 Scalar functions of a matrix === ====== ========= == = ====== int m = A.Nrows(); // number of rows int n = A.Ncols(); // number of columns Real ssq = A.SumSquare(); // sum of squares of elements Real sav = A.SumAbsoluteValue(); // sum of absolute values Real s = A.Sum(); // sum of values Real mav = A.MaximumAbsoluteValue(); // maximum of absolute values Real norm = A.Norm1(); // maximum of sum of absolute values of elements of a column Real norm = A.NormInfinity(); // maximum of sum of absolute values of elements of a row Real t = A.Trace(); // trace LogandSign ld = A.LogDeterminant(); // log of determinant Boolean z = A.IsZero(); // test all elements zero MatrixType mt = A.Type(); // type of matrix Real* s = Store(); // pointer to array of elements int l = Storage(); // length of array of elements 'A.LogDeterminant()' returns a value of type LogandSign. If ld is of type LogAndSign use ld.Value() to get the value of the determinant ld.Sign() to get the sign of the determinant (values 1, 0, -1) ld.LogValue() to get the log of the absolute value. 'A.IsZero()' returns Boolean value TRUE if the matrix 'A' has all elements equal to 0.0. 'MatrixType mt = A.Type()' returns the type of a matrix. Use '(char*)mt' to get a string (UT, LT, Rect, Sym, Diag, Band, UB, LB, Crout, BndLU) showing the type (Vector types are returned as Rect). The versions Sum(A), SumSquare(A), SumAbsoluteValue(A), MaximumAbsoluteValue(A), Trace(A), LogDeterminant(A), Norm1(A), NormInfinity(A) can be used in place of A.Sum(), A.SumSquare(), A.SumAbsoluteValue(), A.MaximumAbsoluteValue(), A.Trace(), A.LogDeterminant(), A.Norm1(), A.NormInfinity(). 3.9 Submatrices === =========== A.SubMatrix(fr,lr,fc,lc) This selects a submatrix from 'A'. The arguments fr,lr,fc,lc are the first row, last row, first column, last column of the submatrix with the numbering beginning at 1. This may be used in any matrix expression or on the left hand side of =, << or Inject. Inject does not check no information loss. You can also use the construction Real c; .... A.SubMatrix(fr,lr,fc,lc) = c; to set a submatrix equal to a constant. The follwing are variants of SubMatrix: A.SymSubMatrix(f,l) // This assumes fr=fc and lr=lc. A.Rows(f,l) // select rows A.Row(f) // select single row A.Columns(f,l) // select columns A.Column(f) // select single column In each case f and l mean the first and last row or column to be selected (starting at 1). If SubMatrix or its variant occurs on the right hand side of an = or << or within an expression its type is as follows A.Submatrix(fr,lr,fc,lc): If A is RowVector or ColumnVector then same type otherwise type Matrix A.SymSubMatrix(f,l): Same type as A A.Rows(f,l): Type Matrix A.Row(f): Type RowVector A.Columns(f,l): Type Matrix A.Column(f): Type ColumnVector If SubMatrix or its variant appears on the left hand side of = or << , think of its type being Matrix. Thus 'L.Row(1)' where 'L' is LowerTriangularMatrix expects 'L.Ncols()' elements even though it will use only one of them. If you are using = the program will check for no loss of data. If you are are using the submatrix facility to build a matrix from a small number of components, consider instead using the concatenation operators [3.6]. 3.10 Change dimensions ==== ====== ========== The following operations change the dimensions of a matrix. The values of the elements are lost. A.ReDimension(nrows,ncols); // for type Matrix or nricMatrix A.ReDimension(n); // for all other types, except Band A.ReDimension(n,lower,upper); // for BandMatrix A.ReDimension(n,lower); // for LowerBandMatrix A.ReDimension(n,upper); // for UpperBandMatrix A.ReDimension(n,lower); // for SymmetricBandMatrix Use 'A.CleanUp()' to set the dimensions of 'A' to zero and release all the heap memory. Remember that 'ReDimension' destroys values. If you want to 'ReDimension', but keep the values in the bit that is left use something like ColumnVector V(100); ... // load values V = V.Rows(1,50); // to get first 50 values. If you want to extend a matrix or vector use something like ColumnVector V(50); ... // load values { V.Release(); ColumnVector X=V; V.ReDimension(100); V.Rows(1,50)=X; } // V now length 100 3.11 Change type ==== ====== ==== The following functions interpret the elements of a matrix (stored row by row) to be a vector or matrix of a different type. Actual copying is usually avoided where these occur as part of a more complicated expression. A.AsRow() A.AsColumn() A.AsDiagonal() A.AsMatrix(nrows,ncols) A.AsScalar() The expression 'A.AsScalar()' is used to convert a 1 x 1 matrix to a scalar. 3.12 Multiple matrix solve ==== ======== ====== ===== If 'A' is a square or symmetric matrix use CroutMatrix X = A; // carries out LU decomposition Matrix AP = X.i()*P; Matrix AQ = X.i()*Q; LogAndSign ld = X.LogDeterminant(); rather than Matrix AP = A.i()*P; Matrix AQ = A.i()*Q; LogAndSign ld = A.LogDeterminant(); since each operation will repeat the LU decomposition. If 'A' is a BandMatrix or a SymmetricBandMatrix begin with BandLUMatrix X = A; // carries out LU decomposition A CroutMatrix or a BandLUMatrix can't be manipulated or copied. Use references as an alternative to copying. Alternatively use LinearEquationSolver X = A; This will choose the most appropiate decomposition of 'A'. That is, the band form if 'A' is banded; the Crout decomposition if 'A' is square or symmetric and no decomposition if 'A' is triangular or diagonal. If you want to use the LinearEquationSolver '#include newmatap.h'. 3.13 Memory management ==== ====== ========== The package does not support delayed copy. Several strategies are required to prevent unnecessary matrix copies. Where a matrix is called as a function argument use a constant reference. For example YourFunction(const Matrix& A) rather than YourFunction(Matrix A) Skip the rest of this section on your first reading. ------------------------------------------------------------------------------ Gnu g++ (< 2.6) users please read on; if you are returning matrix values from a function, then you must use the ReturnMatrix construct. ------------------------------------------------------------------------------ A second place where it is desirable to avoid unnecessary copies is when a function is returning a matrix. Matrices can be returned from a function with the return command as you would expect. However these may incur one and possibly two copyings of the matrix. To avoid this use the following instructions. Make your function of type ReturnMatrix . Then precede the return statement with a Release statement (or a ReleaseAndDelete statement if the matrix was created with new). For example ReturnMatrix MakeAMatrix() { Matrix A; ...... A.Release(); return A; } or ReturnMatrix MakeAMatrix() { Matrix* m = new Matrix; ...... m->ReleaseAndDelete(); return *m; } If your compiler objects to this code, replace the return statements with return A.ForReturn(); or return m->ForReturn(); If you are using AT&T C++ you may wish to replace 'return A;' by return '(ReturnMatrix)A;' to avoid a warning message; but this will give a runtime error with Gnu. (You can't please everyone.) ------------------------------------------------------------------------------ Do not forget to make the function of type ReturnMatrix; otherwise you may get incomprehensible run-time errors. ------------------------------------------------------------------------------ You can also use '.Release()' or '->ReleaseAndDelete()' to allow a matrix expression to recycle space. Suppose you call A.Release(); just before 'A' is used just once in an expression. Then the memory used by 'A' is either returned to the system or reused in the expression. In either case, 'A''s memory is destroyed. This procedure can be used to improve efficiency and reduce the use of memory. Use '->ReleaseAndDelete' for matrices created by new if you want to completely delete the matrix after it is accessed. 3.14 Efficiency ==== ========== The package tends to be not very efficient for dealing with matrices with short rows. This is because some administration is required for accessing rows for a variety of types of matrices. To reduce the administration a special multiply routine is used for rectangular matrices in place of the generic one. Where operations can be done without reference to the individual rows (such as adding matrices of the same type) appropriate routines are used. When you are using small matrices (say smaller than 10 x 10) you may find it faster to use rectangular matrices rather than the triangular or symmetric ones. 3.15 Output ==== ====== To print a matrix use an expression like Matrix A; ...... cout << setw(10) << setprecision(5) << A; This will work only with systems that support the AT&T input/output routines including manipulators. You need to #include the file newmat9.h. The present version of this routine is useful only for matrices small enough to fit within a page or screen width. To print several vectors or matrices in columns use a concatenation operator [3.6]: Matrix A, B; ..... cout << setw(10) << setprecision(5) << (A | B); 3.16 Unspecified type ==== =========== ==== Skip this section on your first reading. If you want to work with a matrix of unknown type, say in a function. You can construct a matrix of type 'GenericMatrix'. Eg Matrix A; ..... // put some values in A GenericMatrix GM = A; A GenericMatrix matrix can be used anywhere where a matrix expression can be used and also on the left hand side of an '='. You can pass any type of matrix (excluding the Crout and BandLUMatrix types) to a 'const GenericMatrix&' argument in a function. However most scalar functions including Nrows(), Ncols(), Type() and element access do not work with it. Nor does the ReturnMatrix construct. See also the paragraph on LinearEquationSolver [3.12]. An alternative and less flexible approach is to use Basematrix or GeneralMatrix. Suppose you wish to write a function which accesses a matrix of unknown type including expressions (eg 'A*B'). Then use a layout similar to the following: void YourFunction(BaseMatrix& X) { GeneralMatrix* gm = X.Evaluate(); // evaluate an expression // if necessary ........ // operations on *gm gm->tDelete(); // delete *gm if a temporary } See, as an example, the definitions of 'operator<<' in newmat9.cxx. Under certain circumstances; particularly where 'X' is to be used just once in an expression you can leave out the 'Evaluate()' statement and the corresponding 'tDelete()'. Just use 'X' in the expression. If you know YourFunction will never have to handle a formula as its argument you could also use void YourFunction(const GeneralMatrix& X) { ........ // operations on X } Do not try to construct a GeneralMatrix or BaseMatrix. 3.17 Cholesky decomposition ==== ======== ============= Suppose S is symmetric and positive definite. Then there exists a unique lower triangular matrix L such that L * L.t() = S. To calculate this use SymmetricMatrix S; ...... LowerTriangularMatrix L = Cholesky(S); If S is a symmetric band matrix then L is a band matrix and an alternative procedure is provied for carrying out the decomposition: SymmetricBandMatrix S; ...... LowerBandMatrix L = Cholesky(S); 3.18 QR decomposition ==== == ============= This is a variant on the usual QR transformation. Start with matrix / 0 0 \ s \ X Y / n s t Our version of the QR decomposition multiplies this matrix by an orthogonal matrix Q to get / U M \ s \ 0 Z / n s t where 'U' is upper triangular (the R of the QR transform). This is good for solving least squares problems: choose b (matrix or row vector) to minimize the sum of the squares of the elements of Y - X*b Then choose 'b = U.i()*M;' The residuals 'Y - X*b' are in 'Z'. This is the usual QR transformation applied to the matrix 'X' with the square zero matrix attached concatenated on top of it. It gives the same triangular matrix as the QR transform applied directly to 'X' and generally seems to work in the same way as the usual QR transform. However it fits into the matrix package better and also gives us the residuals directly. It turns out to be essentially a modified Gram-Schmidt decomposition. Two routines are provided: QRZ(X, U); replaces 'X' by orthogonal columns and forms 'U'. QRZ(X, Y, M); uses 'X' from the first routine, replaces 'Y' by 'Z' and forms 'M'. The are also two routines 'QRZT(X, L)' and 'QRZT(X, Y, M)' which do the same decomposition on the transposes of all these matrices. QRZT replaces the routines HHDecompose in earlier versions of newmat. HHDecompose is still defined but just calls QRZT. 3.19 Singular value decomposition ==== ======== ===== ============= The singular value decomposition of an m x n matrix 'A' (where m >= n) is a decomposition A = U * D * V.t() where 'U' is m x n with 'U.t() * U' equalling the identity, 'D' is an n x n diagonal matrix and 'V' is an n x n orthogonal matrix. Singular value decompositions are useful for understanding the structure of ill-conditioned matrices, solving least squares problems, and for finding the eigenvalues of 'A.t() * A'. To calculate the singular value decomposition of 'A' (with m >= n) use one of SVD(A, D, U, V); // U (= A is OK) SVD(A, D); SVD(A, D, U); // U (= A is OK) SVD(A, D, U, FALSE); // U (can = A) for workspace only SVD(A, D, U, V, FALSE); // U (can = A) for workspace only The values of 'A' are not changed unless 'A' is also inserted as the third argument. 3.20 Eigenvalue decomposition ==== ========== ============= An eigenvalue decomposition of a symmetric matrix 'A' is a decomposition A = V * D * V.t() where 'V' is an orthogonal matrix and 'D' is a diagonal matrix. Eigenvalue analyses are used in a wide variety of engineering, statistical and other mathematical analyses. The package includes two algorithms: Jacobi and Householder. The first is extremely reliable but much slower than the second. The code is adapted from routines in "Handbook for Automatic Computation, Vol II, Linear Algebra" by Wilkinson and Reinsch, published by Springer Verlag. Jacobi(A,D,S,V); // A, S symmetric; S is workspace, // S = A is OK; V is a matrix Jacobi(A,D); // A symmetric Jacobi(A,D,S); // A, S symmetric; S is workspace, // S = A is OK Jacobi(A,D,V); // A symmetric; V is a matrix EigenValues(A,D); // A symmetric EigenValues(A,D,S); // A, S symmetric; S is for back // transforming, S = A is OK EigenValues(A,D,V); // A symmetric; V is a matrix The values of 'A' are not changed unless 'A' is also inserted as the third argument. If you need eigenvectors use one of the forms with matrix 'V'. The eigenvectors are returned as the columns of 'V'. 3.21 Sorting ==== ======= To sort the values in a matrix or vector, 'A', (in general this operation makes sense only for vectors and diagonal matrices) use SortAscending(A); or SortDescending(A); I use the Shell-sort algorithm. This is a medium speed algorithm, you might want to replace it with something faster if speed is critical and your matrices are large. 3.22 Fast Fourier transform ==== ==== ======= ========= FFT(X, Y, F, G); // F=X and G=Y are OK where X, Y, F, G are column vectors. X and Y are the real and imaginary input vectors; F and G are the real and imaginary output vectors. The lengths of X and Y must be equal and should be the product of numbers less than about 10 for fast execution. The formula is n-1 h[k] = SUM z[j] exp (-2 pi i jk/n) j=0 where z[j] is stored complex and stored in X(j+1) and Y(j+1). Likewise h[k] is complex and stored in F(k+1) and G(k+1). The fast Fourier algorithm takes order n log(n) operations (for "good" values of n) rather than n**2 that straight evaluation would take. I use the method of Carl de Boor (1980), "Siam J Sci Stat Comput", pp 173-8. The sines and cosines are calculated explicitly. This gives better accuracy, at an expense of being a little slower than is otherwise possible. Related functions FFTI(F, G, X, Y); // X=F and Y=G are OK RealFFT(X, F, G); RealFFTI(F, G, X); 'FFTI' is the inverse trasform for 'FFT'. 'RealFFT' is for the case when the input vector is real, that is Y = 0. I assume the length of X, denoted by n, is even. The program sets the lengths of F and G to n/2 + 1. 'RealFFTI' is the inverse of 'RealFFT'. 3.23 Interface to Numerical Recipes in C ==== ========= == ========= ======= == = This package can be used with the vectors and matrices defined in "Numerical Recipes in C". You need to edit the routines in Numerical Recipes so that the elements are of the same type as used in this package. Eg replace float by double, vector by dvector and matrix by dmatrix, etc. You will also need to edit the function definitions to use the version acceptable to your compiler. Then enclose the code from Numerical Recipes in extern "C" { ... }. You will also need to include the matrix and vector utility routines. Then any vector in Numerical Recipes with subscripts starting from 1 in a function call can be accessed by a RowVector, ColumnVector or DiagonalMatrix in the present package. Similarly any matrix with subscripts starting from 1 can be accessed by an nricMatrix in the present package. The class nricMatrix is derived from Matrix and can be used in place of Matrix. In each case, if you wish to refer to a RowVector, ColumnVector, DiagonalMatrix or nricMatrix X in an function from Numerical Recipes, use X.nric() in the function call. Numerical Recipes cannot change the dimensions of a matrix or vector. So matrices or vectors must be correctly dimensioned before a Numerical Recipes routine is called. For example SymmetricMatrix B(44); ..... // load values into B nricMatrix BX = B; // copy values to an nricMatrix DiagonalMatrix D(44); // Matrices for output nricMatrix V(44,44); // correctly dimensioned int nrot; jacobi(BX.nric(),44,D.nric(),V.nric(),&nrot); // jacobi from NRIC cout << D; // print eigenvalues 3.24 Exceptions ==== ========== This package includes a partial implementation of exceptions for compilers which do not support C++ exceptions. I used Carlos Vidal's article in the September 1992 "C Users Journal" as a starting point. See customising [2.2] for selecting the correct exception option. See the section on error messages [4] for some notes on the messages generated by the exceptions. The rest of this section describes my exception simulation package. But see the end of the section for controlling the action after an exception has been generated. Newmat does a partial clean up of memory following throwing an exception - see the next section. However, the present version will leave a little heap memory unrecovered under some circumstances. I would not expect this to be a major problem, but it is something that needs to be sorted out. The functions/macros I define are Try, Throw, Catch, CatchAll and CatchAndThrow. Try, Throw, Catch and CatchAll correspond to try, throw, catch and catch(...) in the C++ standard. A list of Catch clauses must be terminated by either CatchAll or CatchAndThrow but not both. Throw takes an Exception as an argument or takes no argument (for passing on an exception). I do not have a version of Throw for specifying which exceptions a function might throw. Catch takes an exception class name as an argument; CatchAll and CatchAndThrow don't have any arguments. Try, Catch and CatchAll must be followed by blocks enclosed in curly brackets. I have added another macro Bounce to mean a rethrow, Throw(). This was necessary to enable the package to be compatible with both my exception package and C++ exceptions. All Exceptions must be derived from a class, Exception, defined in newmat and can contain only static variables. See the examples in newmat if you want to define additional exceptions. I have defined 5 clases of exceptions for users (there are others but I suggest you stick to these ones): SpaceException Insufficient space on the heap ProgramException Errors such as out of range index or incompatible matrix types or dimensions ConvergenceException Iterative process does not converge DataException Errors such as attempting to invert a singular matrix InternalException Probably a programming error in newmat For each of these exception classes, I have defined a member function 'void SetAction(int)'. If you call 'SetAction(1)', and a corresponding exception occurs, you will get an error message. If there is a Catch clause for that exception, execution will be passed to that clause, otherwise the program will exit. If you call 'SetAction(0)' you will get the same response, except that there will be no error message. If you call 'SetAction(-1)', you will get the error message but the program will always exit. 3.25 Cleanup after an exception ==== ======= ===== == ========= The exception mechanisms in newmat are based on the C functions setjmp and longjmp. These functions do not call destructors so can lead to garbage being left on the heap. (I refer to memory allocated by "new" as heap memory). For example, when you call Matrix A(20,30); a small amount of space is used on the stack containing the row and column dimensions of the matrix and 600 doubles are allocated on the heap for the actual values of the matrix. At the end of the block in which A is declared, the destructor for A is called and the 600 doubles are freed. The locations on the stack are freed as part of the normal operations of the stack. If you leave the block using a longjmp command those 600 doubles will not be freed and will occupy space until the program terminates. To overcome this problem newmat keeps a list of all the currently declared matrices and its exception mechanism will return heap memory when you do a Throw and Catch. However it will not return heap memory from objects from other packages. If you want the mechanism to work with another class you will have to do four things: 1: derive your class from class Janitor defined in except.h; 2: define a function 'void CleanUp()' in that class to return all heap memory; 3: include the following lines in the class definition public: void* operator new(size_t size) { do_not_link=TRUE; void* t = ::operator new(size); return t; } void operator delete(void* t) { ::operator delete(t); } 4: be sure to include a copy constructor in you class definition, that is, something like X(const X&); Note that the function 'CleanUp()' does somewhat the same duties as the destructor. However 'CleanUp()' has to do the "cleaning" for the class you are working with and also the classes it is derived from. So it will often be wrong to use exactly the same code for both 'CleanUp()' and the destructor or to define your destructor as a call to 'CleanUp()'. 3.26 Non-linear applications ==== ========== ============ Files solution.h, solution.cpp contain a class for solving for x in y = f(x) where x is a one-dimensional continuous monotonic function. This is not a matrix thing at all but is included because it is a useful thing and because it is a simpler version of the technique used in the non-linear least squares. Files newmatnl.h, newmatnl.cpp contain a series of classes for non-linear least squares, to be extended to maximum likelihood in a later release. Documentation for both of these is in the definition files. Simple examples are in sl_ex.cpp and nl_ex.cpp. These packages do not work with the AT&T [2.3.1] and HPUX [2.3.4] compilers and newmatnl is questionable under Sun 4.0.1 [2.3.6]. 4 Error messages = ===== ======== Most error messages are self-explanatory. The message gives the size of the matrices involved. Matrix types are referred to by the following codes: Matrix or vector 1 Symmetric matrix 5 Band matrix 9 Symmetric band matrix 13 Lower triangular matrix 17 Lower triangular band matrix 25 Upper triangular matrix 33 Upper triangular band matrix 41 Diagonal matrix 63 Crout matrix (LU matrix) 65 Band LU matrix 73 Other codes should not occur. See the section on exceptions [3.24] for more details on the structure of the exception classes. I have defined a class Tracer that is intended to help locate the place where an error has occurred. At the beginning of a function I suggest you include a statement like Tracer tr("name"); where name is the name of the function. This name will be printed as part of the error message, if an exception occurs in that function, or in a function called from that function. You can change the name as you proceed through a function with the ReName function tr.ReName("new name"); if, for example, you want to track progress through the function. 5 Bugs = ==== * Small memory leaks may occur when an exception is thrown and caught. * My exception scheme is not properly linked in with the standard library exceptions. In particular, my scheme may fail to catch out-of-memory exceptions. 6 List of files = ==== == ===== README readme file NEWMATA TXT documentation file NEWMATB TXT notes on the package design KBRIGGS TXT Keith Briggs' list of matrix libraries BOOLEAN H boolean class definition CONTROLW H control word definition file INCLUDE H details of include files and options MYEXCEPT H general exception handler definitions NEWMAT H main matrix class definition file NEWMATAP H applications definition file NEWMATIO H input/output definition file NEWMATNL H non-linear optimisation definition file NEWMATRC H row/column functions definition files NEWMATRM H rectangular matrix access definition files PRECISIO H numerical precision constants SOLUTION H one dimensional solve definition file BANDMAT CPP band matrix routines CHOLESKY CPP Cholesky decomposition EVALUE CPP eigenvalues and eigenvector calculation FFT CPP fast Fourier transform HHOLDER CPP QR routines JACOBI CPP eigenvalues by the Jacobi method MYEXCEPT CPP general error and exception handler NEWMAT1 CPP type manipulation routines NEWMAT2 CPP row and column manipulation functions NEWMAT3 CPP row and column access functions NEWMAT4 CPP constructors, redimension, utilities NEWMAT5 CPP transpose, evaluate, matrix functions NEWMAT6 CPP operators, element access NEWMAT7 CPP invert, solve, binary operations NEWMAT8 CPP LU decomposition, scalar functions NEWMAT9 CPP output routines NEWMATEX CPP matrix exception handler NEWMATNL CPP non-linear optimisation NEWMATRM CPP rectangular matrix access functions SORT CPP sorting functions SOLUTION CPP one dimensional solve SUBMAT CPP submatrix functions SVD CPP singular value decomposition EXAMPLE CPP example of use of package EXAMPLE TXT output from example GNU MAK general make file for Gnu G++ CC MAK general make file for AT&T, Sun and HPUX MS_NT MAK general make file for Microsoft C 8.0 (windows NT) WATCOM MAK general make file for Watcom 10 WATCO_NT MAK general make file for Watcom 10 (windows NT) EX_B MAK make file for example for Borland C 3.1 EX_BC45 MAK make file for example for Borland C 4.5 (windows NT) TMT_B MAK make file for test for Borland C 3.1 TMT_BC45 MAK make file for test for Borland C 4.5 (windows NT) TMT_Z MAK make file for test for Zortech SL_EX CPP example of OneDimSolve routine SL_EX TXT output from example NL_EX CPP example of non-linear least squares NL_EX TXT output from example See the section on testing [2.6] for details of test files. 7 Problem report form = ======= ====== ==== Copy and paste this to your editor; fill it out and email to robert.davies@vuw.ac.nz or Compuserve 72777,656. Version: ............... newmat08 (19 Jan 95) Primary site ........... SimTel Downloaded from: ....... Your email address: .... Today's date: .......... Your machine: .......... Operating system: ...... Compiler & version: .... Describe the problem - attach examples if possible: